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We consider a situation commonly encountered in microfluidics: two streams of miscible liquids 
are brought at a junction to flow side by side within a microchannel, allowing solutes to diffuse from 
one stream to the other and possibly react. We focus on two model problems: (i) the transverse 
transport of a single solute from a stream into the adjacent one, (ii) the transport of the product 
of a diffusion-controlled chemical reaction between solutes originating from the two streams. Our 
description is made general through a non-dimensionalized formulation that incorporates both the 
parabolic Poiseuille velocity proflle along the channel and thermal diffusion in the transverse direc- 
tion. Numerical analysis over a wide range of the streamwise coordinate x reveal different regimes. 
Close to the top and the bottom walls of the microchannel, the extent of the diffusive zone follows 
three distinct power law regimes as x is increased, characterized respectively by the exponents 1/2, 
1/3 and 1/2. Simple analytical arguments are proposed to account for these results. 



I. INTRODUCTION 

Microfluidics is a very promising format to measure 
the dynamics of diflerent processes: diffusion of a so- 
lute [H, 0, 0| , protein folding 0, [1 ,01, k inetics of chem- 
ical reactions 0, S|. ■ . (see Refs. |9l. lid. [lll | for reviews 
on microfluidics). In some of these approaches, one fol- 
lows the reaction-diffusion dynamics of solutes between 
parallel streams. Figure [T] displays the simplest case 
of a Y-junction commonly encountered in microfluidic 
experiments, where two liquids are injected in a main 
microchannel. Since flows are laminar given the small 




FIG. 1: Geometry of the problem. The origin of the axes 
is set at the junction between the two inlets and Z ranges 
between ±d/2. A and B indicates the two arms in which 
the miscible liquids are injected. The liquids may react when 
they interdiffuse into each other. 



length scales and velocities involved, the liquids only mix 
by transverse diffusion, and one can relate the distance 
downstream the channel, to the time elapsed since the 
two streams where put into contact. Essentially, if the 
average velocity in the channel is V, one expects that 
the situation at a distance X from the junction corre- 



sponds to the outcome of diffusion-reaction over a time 
t = X/V, so that for example, simple interdiffusion leads 
to a diffusion zone of width Y ~ {2DX/V)^/^. 

However, in confined geometries the situation is actu- 
ally more complex as pressure-driven flows are heteroge- 
nous (parabolic velocity proflle), resulting in residence 
times that vary with the distance to the walls. Therefore, 
in order to extract from these experiments physical pa- 
rameters such as diffusion coefficients or rate constants, 
it is essential to understand correctly the mass transport 
phenomena in microchannels. 

In this context, Ismagilov et al. have studied using 
confocal microscopy, the formation and the diffusion of a 
fluorescent product in the diffusion cone of two reactive 
miscible solutions flowing side by side in a microchannel 
(the geometry is the same as displayed in Fig. [1]) . They 
have shown, both theoretically and experimentally, that 
the extent of the transverse diffusive zone near the top 
and the bottom walls {Z = ±d/2) scales as F ^ X^/^, 
i.e. follows a 1/3 power law in the streamwise coordi- 
nate X. In the central plane of the microchannel, they 
observed the classical 1/2 power law for diffusive pro- 
cesses Y ^ X^l"^ . Such behaviours are due to the cou- 
pling of transverse diffusion and velocity gradients in the 
channel. These measured exponents can be accounted 
for by arguments derived from the classical Leveque prob- 
lem [l3|- Kamholz et al. have performed numerical simu- 
lations of the same problem [1J|, and found the different 
regimes expected. However, their results concern only 
speciflc sets of dimensionalized parameters which makes 
it difficult to apply them to other experimental config- 
urations. In addition, some of the results displayed in 
Ref. are inconsistent with findings shown later in 
the present manuscript. More recently, a significant im- 
provement has been brought by Jimenez who performed 
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numerical simulations of the same problem using dimen- 
sionless equations His work accounts for the differ- 
ent regimes observed experimentally, and the author also 
clarifies these results using analytical solutions of the con- 
centration profile in the range of small and large X. In 
a recent paper (Supporting Information of Ref. [§]), we 
have studied the effect of the Poiseuille flow developped 
in a microchannel of high aspect ratio, on the reaction- 
diffusion dynamics occurring in the diffusion cone of two 
reactive solutions injected side by side. Through numeri- 
cal simulations, we showed that for the high aspect ratio 
studied, a 2D description in terms of height-averaged con- 
centrations and velocities was appropriate and could be 
safely used to extract the rate constant of the reaction 
from experimental data. 

Our aim in this paper is to revisit formally these trans- 
verse transport problems so as to provide generally ap- 
plicable insights as to the outcome of experiments per- 
formed in arbitrary conditions (flow rate, geometry,. . . ), 
but for the requirement that the Peclet number be large 
(see below). In the first part of this paper, we discuss 
the simple case of the transverse diffusion of a solute, us- 
ing numerics and simple analytical arguments in a generic 
non-dimensionalized formulation. We clarify that the 1/3 
power law behaviour for transverse spreading {Y ^ X^^^) 
holds after the junction in layers close to the walls, the 
thickness of which increases as ~ X^^^. For a small but 
finite distance from the walls (i.e. for a given Z), the 
transverse extent Y of the diffusion zone actually follows 
successively three power law regimes in X as X is in- 
creased, with exponents 1/2, then 1/3 and finally 1/2 
again. We then discuss the relevance of these findings 
for published data. In a second part, we study the more 
complex case of coupled diffusion and reaction of two re- 
active solutes yielding the apparition and diffusion of a 
product in the neighbourhood of the interface between 
the two streams. Most of the phenomenology described 
in the first part is shown to persist when the kinetics of 
the reaction is fast compared to diffusion, with however 
a few qualitative differences that are pointed out. 



II. DIFFUSION OF A SOLUTE INTO A 
NEIGHBOUR STREAM 

A. Assumptions of the model 

We consider a microchannel with a high aspect ratio 
r = L/d^ 1, as depicted in Fig. [TJ Two solutions of 
the same solvent are injected at a same constant fiow 
rate in the two arms of the microchannel: (i) the solu- 
tion fiown through arm A contains at a dilute level a 
solute A, whereas (ii) pure solvent is fiown through the 
other arm B. We are interested in the concentration pro- 
files A{X, Y, Z) of this solute downstream. Given our mi- 
crofiuidic motivation, we consider that the fiow is strictly 
laminar (small Reynolds number) and that the solute 
only disperses by molecular diffusion. We also suppose 



that the diffusion coefficient D of the solute, as well as 
the fiuid viscosity and density, are constant for the dilute 
solutions considered here. 

With these assumptions, the velocity profile is locally 
parabolic and well described by Hele-Shaw formulas (ex- 
cept for positions Y at a distance d or smaller from the 
two lateral walls at Y — ±L/2). Typically after an en- 
trance region of order L from the junction, the fiow takes 
the simple form v{X, Y, Z) = u(Z)X with 

with U the maximal velocity. 

We now write the equation describing the transport of 
the solute resulting from convection in this velocity field 
and thermal diffusion. We focus immediately on the limit 
of high Peclet numbers, i.e. Pe = J7c?/Z) » 1, where one 
can neglect the diffusion along the fiow direction, and the 
concentration A of solute A evolves according to 

v{Z)dxA{X, r, Z) = D{dl + dl)A{X, Y, Z) . (2) 

If we neglect entrance effects, i.e. if we assume that the 
solute does not significantly diffuse transversely in the en- 
trance zone, the concentration downstream can be found 
by solving Eq. |[2]) with no-fiux boundary conditions on 
the lateral walls, and "initial" conditions at X ~ 0: 
A{Y, Z) = for y > and A{Y, Z) = A^ for Y <Q. 
We make the problem non-dimensional using 

X = X/(dPe), z = Z/d, y = Y/d, and a = A/A^, (3) 

so that Eq. ^ takes a parameter-free form 

(1 - {2zf) d.,a{x, y, z) = (9^ + dl)a{x, y, z) . (4) 

B. Numerical computation of the model 

We solve Eq. ^ using a simple numerical approach 
(the Euler method). We descretize the (y, z) plane using 
a 2D "grid" with nodes yj = -T/2 + dy/2 + (j - l)dy, 
and Zj — —1/2 -I- Az/2 + [j — l)dz, with Ay = T/n and 
dz = \/p, with typically p = 60 points across the channel 
height and n = 80 points across the channel width. We 
impose the classical no-fiux boundary conditions at the 
walls at z = ±1/2 and y = ±T/2. As in Ref. [l3|, we 
perform different simulations and overlap their outcome 
to cover a wide range of values of x. More precisely, 
to reach sufficient accuracy over this whole range, the 
resolution on the y-axis is varied through variations of 
r, each simulations being terminated before the solute 
diffuses up to the side walls as y = ±r/2. 

C. Numerical results 

Typical concentration profiles within cross-sections at 
different downstream locations x are shown in Fig. [2l 



FIG. 2: Typical z-y slices of the concentration a of the solute 
at (a) X = 10"*^, (b) 10~^, and (c) 10~^. A linear gray scale 
is used to code the values of a: black corresponds to 0, and 
white corresponds to 1. 

For X < 0.1, the concentration profiles are not homoge- 
nized across the height of the microchannel, as already 
explained in Refs. [13, . This is a consequence of 
the dispersion of the residence times in the microchannel 
due to the Poiseuille fiow. For x > 0.1, one almost recov- 
ers homogeneous profiles. This is expected since the cri- 
terion for the solute to sample all z-positions by diffusion 
reads roughly a; > 1/8, as the solute has to diffuse over 
the half height of the channel so that 2DX/V > {d/2f. 

To get further insight into the broadening with x of the 
concentration front, we define an average width 'w{x, z) 
of the concentration profiles, at height z and distance x, 
as 

nT/2 .r/2 
2w'^{x,z) = / dyy^dya{x,y,z)/ / dy dya{x,y, z) . 
J-r/2 J-r/2 

(5) 

In the ideal case of simple diffusion with a ho- 
mogeneous velocity profile (plug fiow), a{x,y,z) = 
0.5 (1 -I- erf(— y/(2ii;(x))), and the diffuse broadening 
obeys the simple diffusive scaling w{x) ~ x^/^ with an 
exponent 1/2. 

The curves z) vs. x for a Poiseuille fiow profile are 
shown in Fig. [3] for several heights z in the microchan- 
nel. For values of x below 0.1, a dispersion of the widths 
w{x, z) is visible, as expected from Fig. [51 the concentra- 
tion profiles are wider close to the bottom and the top 
walls. For values of x greater than 0.1, all the curves 
collapse and a homogenous (along z) diffusion along y is 
recovered. 

To describe in more detail the /oca/ behaviour of w(a;, z) 
vs. x, we fit these curves locally by power laws w{x, z) ~ 



FIG. 3: lo 

Sio(™) loSio(^') for different z. The arrow 
indicates the range of increasing z from —l/2-\-6.z/2 to —Az/2 
{Az = 1/60 in the present simulations). The two continous 
lines indicate the 1/2 and 1/3 power law behaviours. 

^oi[x,z) ^ sliding decade of x, from x ~ 10~^ to x « 1. 
The exponents a{x, z) found by such an analysis are dis- 
played in Fig.m For very low values of a; < 10^'', one ob- 




FIG. 4: a{z,x) vs. logj^Q(a;) and z. The exponents a{x,z) 
are estimated by local fits of the curves displayed in Fig. [3] 
by power laws. The dashed white line indicates u — 1 [see 
Eqs. l[6ll- lll0p ]. and the continous white line corresponds to 
Eq. lf2T|l with Pe — 1000. For clarity, these two lines have 
been only plotted for z > 0. 

serves that w{x, z) ~ x^/^ for all the investigated range 
of z. For 10~^ < a; < 10~^, the curves corresponding 
to z-positions close to the walls reach a 1/3 power law 
regime. In this regime, the local exponent in the middle 
of the microchannel deviates from 1/2, reaching values 
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up to 0.55. Eventually, the classical 1/2 exponent of ho- 
mogenous diffusion is recovered for x > 0.1. 

At small but finite distance from the walls (< 0.2), we 
therefore predict a succession of three power law regimes 
with exponents 1/2, 1/3, and 1/2. The transition from 
the 1/2 to the 1/3 regime occurs earher for values of z 
closer to the walls. To our knowledge, no experiments 
have shown this transition. 

In many experiments, the quantity measured is the 
concentration field averaged over the thickness of the 
microchannel. We compute this quantity from our so- 
lutions of Eq. (|4]), and define the corresponding width 
Wm (a^) describing broadening in this averaged map using 
the analog of Eq. |[5]). We can now quantitatively an- 
swer one of the main points raised in the introduction: 
for very thin channel, how good/bad is the approxima- 
tion consisting in assuming fast diffusion over Z so that 
flow heterogeneities can altogether be neglected ? This 
approximation le ads to the homogeneous diffusion pre- 
diction w{x) ~ ^J?>x/2. We compare it to the real width 
of the averaged map Wm[x) by plotting in Fig. [UJa) a 
measure of the relative err or ma de in using this approx- 
imation 5{x) = Wrn{x) / — 1. From the plot it 
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FIG. 5: (a) 5 — Wm{x)/ ^3x/2 — 1 vs. x. 5 allows one 
to estimate the deviation between the width of the averaged 
concentration profiles and the solution when the Poiseuille 
flow has not been taken into account, (b) SW vs. x. 5W 
is the difference between the widths of the diffusion zone at 
z — > ±1/2 and at 2 = 0. 

is clear that while important deviations are observed at 
very short distances, the error is less than 5% for all val- 
ues of X larger than 0.03, which corresponds in real units 
to distances after the junction X < 0.03 Fed. 

Figure \Mp) shows another measure of the heterogene- 



ity along z of the concentration profile: 5W, the differ- 
ence between the widths of the interdiffusion zone at the 
wall and at z = 0. As can be seen on this plot, the max- 
imal difference is reached for x « 10^^, and corresponds 
to values of about 7% (in real units 0.07 d). 



D. Transition of the exponent from 1/2 to 1/3 
close to the walls 

We present now analytical arguments that help under- 
stand the above results, namely the transition from the 
1/2 to the 1/3 regime observed close to the walls, for a 
value of x that increases with the distance to the wall. 

We focus on the vicinity of the bottom wall, and make 
therefore the additional reasonable approximation that 
the velocity profile there is almost linear, so: v(Z') = 
jZ', where 7 = AU/d is the local shear rate at the wall 
and Z' = Z + d/2 is the distance to the wall. In the dif- 
fusion process, the only relevant length scale is therefore 
I = y^D/j. Close to the wall, the concentration field of 
the solute A evolves through 

Z'dxa=^l^{d^ +dl,)a. (6) 

We make a change of variables akin to the classical treat- 
ment of the Leveque problem: 

u = Z7(Xi/3p/3)^ (7) 

V = y/(Xl/3/2/3) , (8) 

a = Fiu,v), (9) 
so that Eq. ^ reads 

{dl + d'-)F^-^u[uduF + vd,F]. (10) 
The boundary and initial conditions 

{dz'a)z'=o^O, (11) 
a = for r > 0, and a = 1 for Y <0, (12) 

read now 

{duF)u=o = , (13) 
lim F = 1 , (14) 

lim F = Q. (15) 

V — V — 00 

Using Eq. (HI, the local width of the concentration profile 
is: 

w\X,Z')=X^/H''^'^G{u), (16) 

where 

G{u) = / dvv^d.F. (17) 
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Using Eq. (fTO|) . G{u) has to satisfy 
-2 2u 



G" 



-G' 



G + 2 = 0. 



with 



lim G . 



(18) 



(19) 
(20) 



Close to the wall Z' ~ and u ~ 0, the conditions 
{duG)u=o = therefore gives G(u) = cste and w{X, Z) ~ 
X^/^Z'^/'^. This is the classical results of the Levfique 
problem [l^ and shown by our simulations for positions 
Z' very close to the walls, for x < 0.1. For large u, 
one should have w'^{X, Z') = KDX/{jZ'), and therefore 
G(m) = K/u. Using Eq. ([18]), one has K = 2 which is 
the classical diffusive behaviour. 

Figure El helps to sum up these results. At a given po- 
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FIG. 6: Schematic representation of the asymptotic behaviour 
of the solutions of Eq. ifTS]! . The line X = corresponds to 
?i ~ oo for finite Z' . The dashed line indicates u = 1. 

sition Z' , one encounters two regimes as X is increased. 
For X ~ 0, one has u -> oo and w{X,Z') - X^/^, 
and for large X one has u — > and 'w{X,Z') ^ X^^^. 
The boundary between these two regimes can be esti- 
mated using the curve u = 1, i.e. Z' = X^^^P^^. In 
the previous systems of units [see Eq. u = 1 reads 
z ~ ±0.5 =F {x/4y^^. One of these two boundaries is 
shown in Fig. [4] and a good agreement is observed with 
the frontier between the 1/2 and 1/3 regimes. 

Note that this simple argument, that permits to appre- 
hend the transition from 1/2 to 1/3, is strictly speaking 
limited to the vicinity of the walls. Further away, one 
should take into account the curvature of the flow pro- 
file, and eventually describe the merger of the two zones 
of 1/3 behaviour. 



E. Validity of the approximations embedded in the 
model 



We have altogether neglected longitudinal diffusion 
(along X) on the ground of a large Peclet number Pe = 



Ud/D ^ 1. To assess more locally the range of validity 
of this approximation, a useful quantity is v{Z)X/D, the 
Peclet number using the length scale X and the velocity 
v{Z). Indeed, positions X,Z such that v{Z)X/D ^ 1 
correspond to the region where diffusion along the ax- 
ial direction is negligible. The curve v{Z)X/D = 1 is 
thus a reasonable estimate for the boundary of the do- 
main where our approximation should be valid. In non- 
dimensionalized units, this boundary reads 



(21) 



Pe^ 1 - (2z)2 



and is plotted in Fig. |4] for the case Pe = 1000. It inter- 
sects the curve u = 1 indicating the growing boundary 
layer corresponding to the 1/3 power law regime for val- 
ues X w 1/Pe^ and z ± 1/2 w ±{l/AVe^f'^. So, to vi- 
suaHze experimentally comfortably this 1/3 regime, one 
should use high Peclet numbers, and focus on positions 
(z,x) given by the previous relations. For Pe = 1000, 
the crossing point (see Fig. H]) is at z ± 1/2 « 610"'^, 
i.e. below 1% of the height of the microchannel, leaving 
ample space for observing this regime. 

Another issue concerns having neglected any specific 
effect in the entrance region of size L. In this zone, the ve- 
locity profile evolves from what it was in each branch into 
its final form given Eq. (1). The extent of this region in 
non-dimensional units is of the order of x < (L/d) Pe~^, 
so that unless Pe is really very large, this region will over- 
lap with the 1/3-regime. For example for Pe = 1000 and 
L/d = 10, the entrance region corresponds to x < 0.01. 
However, we do not think that a refined description with 
a model for this entrance region would modify the pic- 
ture presented. Indeed, it takes only a distance of order 
X ^ d after the apex of the junction to reach a parabolic 
profile in the interfacial region, and it is only the ampli- 
tude of this parabola (the local maximal velocity) that 
relaxes to its asymptotic value over X ^ L. 



F. Comparison with experiments 

We now compare our results with existing experimen- 
tal data. Kamholz et al. have measured the diffusivity of 
various fiuorescent molecules using the device sketched 
in Fig. [T] [l| . They used microchannels of small thickness 
(d = 10 /im) and models assuming fast averaging over 
this thin dimension to extract the diffusion constants. 
Our analysis shows that this procedure was actually ap- 
propriate given its simplicity. Indeed, their "worst" ex- 
perimental configuration (highest fiow rate ~ 1000 nL/s, 
and smallest diffusivity D = 6.2 10"^^ m^/s for strepta- 
vidin, so Pe w 8000), corresponds to measurements taken 
at a; w 0.05 in non-dimensionalized units, for which the 
relative error from using the simplified model is less than 
5% as can be seen on Fig. [D^a) (even if the difference of 
width 5W is maximal at this position [see Fig. [UJb)]). 

As pointed in the introduction of the present pa- 
per, our results are not consistent with the theoreti- 
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cal data discussed in Ref. In this work, the au- 

thors do no use non-dimensionaUzed equations but have 
computed the concentration profiles for a specific case 
(flow rate Q = 42 nL/s, L = 2405 /im, d = 10 /im, 
D = 3.4 10-1° s-\ leading to Pe « 50). They found, 
in the midplane of the channel (z = 0) and in the range 
X = 10-^-10"^, a power law regime for the transverse 
spreading given by an exponent 2/3. Such a result is 
not consistent with the data displayed in Fig. IH and we 
believe it is due to the fitting procedure of the curves 
w{x, z) vs. X. More important, the authors claim that 
a model with a homogeneous velocity profile (ID model) 
overpredicts the diffusion width by ~ 35%, even for dis- 
tances X = 1000 iim, i.e. a: = 2 in non-dimensionalized 
units. Such results are in contradiction with our data 
[see Fig.ISlJa)]. 

We now consider the data of Ismagilov et al. [l^ . 
who used confocal fiuorescence microscopy to measure 
the transverse diffusive broadening of a fiuorescent probe 
formed as the product of a reaction involving diffusing 
reactants originating from solutions fiowing side by side 
in a microchannel. Clearly, that situation is more com- 
plex than the one analyzed above as both diffusion and 
reaction come into play. The data displayed in Ref. [l^ . 
more precisely in their Fig. 2, concern a microchannel 
of thickness d = 105 /im, a maximal velocity of the or- 
der of 12 cm/s, and an investigated range of X ranging 
between 10^ and 10^ fjim. Taking D k, 10^*^ m^/s to es- 
timate Pe ~ 12600, the dimensionless x ranges roughly 
between 10^"* and 10^^. As can be seen on Figs. [3] and HI 
our analysis suggests that in this range, the boundary 
layers where the transverse diffusive zone widens as x^^^ 
should occupy a significant portion of the microchannel, 
in agreement with the observation by the authors of such 
a scaling for the broadening close to the walls. However, 
as can be seen on Fig. [5l[b), from our analysis, the dif- 
ference between diffusive width at the walls and at the 
center should be less than « 0.07c?, i.e. « 7 /im. The 
values reported in Ref. lead to a somewhat larger 
estimate (~ 30 /im)- the difference is too large to be 
simply the result of a slightly different definition of the 
diffusion width w [see Eq. (O]-. Furthermore, the inter- 
diffusion zones they observed have a different shape than 
those displayed in Fig. [2l We conclude that while the 
presence and location of a "1/3 regime" is robust to the 
differences between the two situations, some features do 
show a difference. To confirm this we now proceed to a 
brief investigation of a reaction-diffusion model. 



III. DIFFUSION OF QUICKLY REACTING 
SPECIES 

A. Assumptions of the model 

We turn to the case of a chemical reaction A -I- B ^ C. 
Different groups have theoretically investigated the solu- 
tions of the unidimensional problem of reaction-diffusion 



without advection (l6l. [TtI. fT^ . These works mainly deal 
with irreversible reactions A-l-B — > C, and have provided 
scaling laws of the width of the reaction front, both in 
the asymptotic and in the short-time limits. In our case, 
we are interested in the concentration map of the prod- 
uct C, when the two reactants, A and B, are injected in 
the two arms of the microchannel (see Fig.[l|), at a same 
given constant fiow rate, and at concentrations oq and 
6o respectively. As above, we restrict ourselves to dilute 
solutions, and further focus on fast chemical reactions, 
so that chemical equilibrium AB/C = if eg is satisfied at 
each point of the microchannel {K^q is the equilibrium 
constant of the reaction, and A, B and C are the concen- 
trations of species A, B, and C respectively). To reduce 
the number of independent parameters, we take the diffu- 
sion coefficients of species B and C to be almost identical, 
i.e. Dc ~ Db- This should well describe complexation 
reactions where A is a small molecule that can not affect 
the diffusion coefficient of the larger B. Reaction-diffusion 
dynamics in a microchannel with the same assumptions 
than previously read 

v{Z)dxA =-R + DA{dl+dl)A, (22) 
v{Z)dxB =-R + DB{dl+dl)B, 
v{Z)dxC ^+R + DB{d^+dl)C , 

R corresponds to the reaction term. Since the reaction is 
very fast, AB/C = Keq, and the problem is governed by 
only two equations. Considering species conservation, we 
focus on the combination: P = A + C, and Q = B + C, 
so Eqs. (I22|) read 

viZ)dxP = DA{dl- + dl)P (23) 

+ {DA-DB){d'^+dl)C , 
v{Z)dxQ = DB{d^+dl)Q, 

where C = C{P, Q) is to be understood as the implicit 
solution the local chemical equilibrium: 

KeqC ^ {P - C){Q - C) . (24) 

We move to dimensionless variables, 

p^P/aa, q = Q/ba, c=C/{aabof'\ (25) 

and again : 

x = X/{dVe), z^Z/d, and y^Y/d, (26) 
with Pe = Ud/DA, Eqs. ([23]) and ^ read 

{l-(2zf)d.p=idl+dl)ip+^^c), (27) 

{l-i2zf)d,q^x\dl+d^_)q, 
- c(7 + f3p + q/f3) +pq = 0, 

where we now have three dimensionless parameters that 
control the outcome: 

= Db/Da, (i^ = ao/feo, 7 = K,q/{aoboY'^ . (28) 
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As can be seen on Eqs. I|27p . the case x = 1 yields p = q 
and is equivalent to the problem discussed previously. In 
general X 7^ 1 and /3 7^ 1, so the diffusion pattern is not 
symmetric. 



B. Numerical computation and results 

As above, we write Eqs. ([27l) on a discrete set of n = 
100 X p = 80 points. We impose the classical no-flux 
boundary conditions at all the walls of the microchannel, 
and the initial conditions at X = are p = 1 {q = 
resp.) for Y > 0, andp = {q = 1 resp.) for Y <0. 

Figure [7] shows the results obtained for a speciflc set 
of parameters x ~ 0-3, /? ~ V200, and 7 = 0.0055, 
chosen because it corresponds rou ghly to the situation 
studied by Ismagilov et al. in Ref. |12| {Da ~ 1.2 lO'^, 




FIG. 7: (a) concentration field c of the product C averaged 
over the height of the microchannel in the case of a chemical 
reaction described by Eqs. lf27|l with x ~ 0.3, (3 — \/200, and 
7 — 0.0055. (b-c) two slices y vs. z of the concentration of 
the product of the reaction, at the locations specified by white 
lines on (a). 

Db w 10-1" m^s-i, ao 1 mM, 60 = 5 /iM, K^q = 
0.39 pM). The depth-averaged concentration map c in 
the microchannel is displayed in this figure (fiow is down- 
wards), together with two cross-sectional plots of the 
concentration in the {y.,z) plane. As for the diffusion 
of a solute, the diffusion zone is not homogeneous over 
the height of the channel, due to the presence of the 
Poiseuille profile. Moreover, the diffusion zone is not 
symmetric {x,y), as a consequence of the stoechiometric 



imbalance (/? 7^ 1), and the asymmetry in diffusion coef- 
ficients (x 7^ 1) H- The good agreement between the ex- 
perimental shapes of the concentration profiles obtained 
in Ref. [l3| and our simulation is clear. Moreover, the 
differences 5W between the widths of the interdiffusion 
zone at the wall and at z = 0, are now consistent with 
the experimental data (for example, in Fig. 2 of Ref. [l^ . 
5W w 0.1-0.3, for Pe w 12600 and x w 0.0001-0.001, our 
simulations give roughly the same results). 

As above (see Figs. [3] and H]), we perform numerical 
simulations of Eqs. l(27|) over a wide range of x, and com- 
pute w{x,z). The corresponding results are shown in 
Fig.IH together with a map of the local exponents a{x,z), 
estimated by fitting the curves over a sliding decade of x 
with a law w{x, z) ~ (in this case we have chosen 

n = 80 and p = 40 for the simulations). Obviously the 




FIG. 8: log^Q(?i)) vs. logio(x ) for different z in the case of 
Eqs. ((22l) with x = 0.3, P = \/200, and 7 = 0.0055. The two 
continous lines indicate the 1/2 and 1/3 power law behaviours. 
Insert: corresponding exponents a{x,z) vs. logiQ(a::) and z, 
obtained by local fits of w vs. x by power laws. 

phenomenology is very similar to that for the simple dif- 
fusion of a solute studied earlier: (i) close to the walls at 
z = ±0.5, w ~ x^/'^, whereas for z « 0, then w ~ a;^/^; 

(ii) beyond x ~ 0.1, the diffusion zone is homogeneous 
over the height of the channel and widens as u; ~ a;^/^; 

(iii) close to the walls the sequence of power law regimes 
with exponents 1/2, 1/3, and 1/2 is recovered; (iv) the 
locus of the transition from the 1/2 to the 1/3 regime 
occurs earlier for z closer to the walls. We find similar 
results when solving numerically Eqs. (|27| for different 
set of parameters x, P and 7. However, the increased 
complexity of the equations does not permit to perform 
as easily as in subsection II-D a simple analytical argu- 
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ment to formally argue for this phenomenology to hold 
whatever the values of the parameters. Theoretical ap- 
proaches such as the ones described in Refs. [H, [3| 
may be useful to answer these questions. 

IV. CONCLUSION 

The present numerical and theoretical study shows the 
generality of features reported for specific conditions by 
previous authors on lateral transport of solutes between 
streams in pressure-driven flows. For diffusion alone, the 
variations of the velocity along the thickness of the chan- 
nel lead to remarkable features: while in the midplane 
of the microchannel the diffusive width classically grows 
with a 1/2 exponent, close to the wall the broadening 
downstream can go through three power law regimes with 



exponents 1/2, 1/3 and 1/2 again. Remarkably, our nu- 
merical computations suggest that this phenomenology 
may also hold for reaction-diffusion problems, at least in 
the diffusion-controlled limit of fast reaction kinetics. All 
our discussion and arguments, together with the various 
plots presented, are proposed in non-dimensional form, 
which hopefully should help experimentahsts in assess- 
ing a priori the relevance of 3D effects in microfluidic 
experiments involving the lateral transport of solutes. 
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